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ABSTRACT 

We present two-dimensional hydrodynamical simulations of slowly rotating gas that is 
under the influence of the gravity of a super massive black hole and is iiTadiated by a thin UV 
accretion disc and a spherical X-ray corona. We calculate the accretion luminosity of a system 
based on the accretion-rate which is assumed to be equal to the mass-supply rate at the radius 
of ~ 10~^ pc. For the models with high temperature gas at large radii (~ 10 pc) and high 
luminosities, we find a strong correlation between the mass-outflow rate (A/out) and the lumi- 
nosity (L). The power law index {q) describing the i\/out--^ relation is q = 2.0 (±0.1), which 
is very similar to that for radiation-driven stellar and disc wind models. More surprisingly, 
for high density at large radii, we find steady state solutions with the accretion luminosity ex- 
ceeding the Eddington limit. The super-Eddington accretion proceeds in the equatorial region 
and is possible because the radiation flux from the disc is significantly reduced in the equato- 
rial direction due to the geometrical foreshortening effect. In all models, an outflow is driven 
from an inflow with sub-Keplerian rotation. For high temperature at large radii, the inflow 
occurs over a wide range of the polar angles, and the outflow occurs in a relatively naiTow 
polar cone. However, for the super-Eddington cases with low temperature at large radii, the 
inflow persists only very close to the equatorial plane, resembling a thin accretion disc, while 
the outflow arises in a wide range of radii and polar angles. The geometry of this extreme 
inflow-outflow solution is very similar to a radiation-driven wind from a luminous Keplerian 
accretion disc. For the cold super-Eddington solutions, 7\/out is only very weakly correlated 
with L, i.e. ^ q ^ 0.2. This weaker correlation is mainly caused by a mismatch between 
the direction of escaping photons and the inflowing gas: the radiation is emitted mostly in the 
polar directions whereas the inflowing gas occurs mainly in the equatorial region. 
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1 INTRODUCTION 

Induced by accretion of gas on to a super massive (10^-10^*^ 
black hole (SMBH), active galactic nuclei (AGN) release a large 
amount of energy by radiating photons (10^''-10^''Lq). These 
photons influence the physical properties (e.g. the ionization struc- 
ture, gas dynamics and density distribution) of the vicinity of AGN, 
the AGN host galaxies, and the inter-galactic media of galaxy 
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portance of the feedback from AGN has been recognised in 
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cally important (e.g. Iciotti & Ostrikej 1 1991 1200 ll. |2007| : 
ICiotti et al.l l2009h . The radiation can produce mechanical 
feedback via ra diation pressure on matter in the form of out- 
flows/winds (e.g. Pro gal 20071: Proga. Ostriker. & Kurosavy3 20081: 



Dorodnit syn. Kallman. & Progal l2008al lbl: Kurosavya & Progal 
2008., .20091) . A strong support for the existence of such 
(radiation-driven) winds is often found in the highly 
blueshifted broad absorption line feat ures seen in the 
UV and optical spectr a of A G N (e.^. [ Murray et al] 1 19951; 
Proga, Stone, & Kallmaj I2OO0I: IChartas, Brandt, & Gallagher 
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In our previous st udies l lProgj I2OO7I : [Proga et all [20081 : 

[Kurosawa & Proga[ [2009I . hereafter Papers 1, 11 and 111, respec- 
tively), we performed hydrodynamical simulations of radiatively 
driven AGN outflows produced from the infalling gas. Our mod- 
els include radiation forces due to electron scattering and line pro- 
cesses, and radiative cooling/heating. These models use a fixed ac- 
cretion luminosity during the whole simulations regardless of the 
mass inflow rate through the model inner boundary. In principle, 
the rate at which the gas is supplied to the region below the inner 
boundary (~ 10^^ pc) would influence the mass accretion rate on 
to the SMBH at the centre and its luminosity. This could results 
in models with an inconsistent amount of the radiative feedback 
for a given amount of the mass supplied rate to the inner region 
of AGN. Here we relax this assumption of the constant luminos- 
ity, and instead dynamically couple the accretion luminosity to the 
mass supply rate. With this coupling, we examine the following two 
key issues in this study. 

Firstly, we investigate how the mass-outflow rate depends on 
the self-consistently determined accretion luminosity (or equiva- 
lently the Eddington ratio F of the system). This provides a useful 
information on the mass supply rate and the mass inflow morphol- 
ogy which can be incorpo rated with a smaller scale disc accretion 
model (e.g. [Ohsugij[2007[) . Similarly, the mass-outflow rates pre- 
dicted here can be incorporate d with a larger scale galaxy evolution 
model (e.g. [Ciottietal][2009l) . 

Secondly, we test whether steady state solutions with a super- 
Eddington luminosity are possible when the main radiation source 
has a disc geometry instead of a spherical geometry. The system 
with the disc radiation field has a possible channel of a high mass 
supply rate along a equatorial plane at a large distance from the 
disc itself. Because of the disc geometry, the radiation would peak 
strongly along the pole or disc axis directions, and it would gradu- 
ally decrease toward equatorial direction hence the radiation force 
is much weaker in the equatorial plane. 

In the following section, we describe our method and model 
assumptions. We give the results of our hydrodynamical simula- 
tions in Section [3] Discussions on accretion flows with a super- 
Eddington luminosity and the effect of artificially limiting lumi- 



nosity are given in Section|4] Finally, the conclusions of this study 
are summarised in Section|5] 



2 METHOD 
2.1 Overview 

We m ainly follow the m etho ds used for the ax i symm etric mod- 
els bv lProga etZI j2000l) and [Proga & KallmanI l l2004h . but with 
a change in the treatment of the accretion luminosity (Section r2.4b . 
The model geometry and the assumptions of the SMBH and the 
disc are very similar to those in Papers 1, 11 and 111. Readers 
are referred to Fig. 1 of Paper 111 for the diagram of our ba- 
sic model configuration. A SMBH with its mass Mbh and its 
Schwarzschild radius Rs = 2GA/bh /<? is placed at the centre of 
the polar coordinate system (r, 9). The X-ray emitting corona re- 
gions is defined as a sphere with its radius r, which surrounds the 
SM BH. The geometrically thin and optically thick accretion disc 
(e.g. [Shakura & Sunvaev[|l973h is placed on the equatorial plane 
(6 — n/2 plane). The normal vector of the disc coincides with the 
symmetry axis. The hydrodynamic simulations will be performed 
in the polar coordinate system (assuming axisymmetry) with r be- 
tween the inner boundary ri and the outer boundary r-Q. The ra- 
diation forces, from the corona region (the sphere with its radius 
r*) and the accretion disc, acting on the gas are assumed to be ra- 
dial only. The point-source like approximation for the disc radiation 
pressure is used here since the accretion disc radius (ru) is assumed 
to be much smaller than the inner radius, i.e. r-D <S ri. 

Because of the assumed flat accretion disc geometry, the disc 
radiation flux (jT^isc) peaks in the direction of the disc rotational 
axis, and it gradually decreases as the polar angle increases, 
i.e. JTdisc oc I cosO\. The flow is also irradiated by the corona which 
is assumed to be spherical. Further, we assume that the total accre- 
tion luminosity L consists of two components: (1) Ldisc = /disci 
due to the accretion disc and (2) L* = /,L due to the corona 
(i.e. /disc + f* ~ !)• We assume that the disc emits only UV pho- 
tons, whereas the corona emits only X-rays, i.e. the system UV lu- 
minosity, I/uv = fvvL = Ldisc and the system X-ray luminosity, 

= fxL = L, (in other words /uv = /disc and /x = ./*). 

We include the effect of radiation forces due to electron scat- 
tering and line processes (line force). To evaluate the line force, we 
follow the method described bv lProga et al.l j2000l) who use a modi- 
fied ve rsion of the formalism developed bv^ Caston Abbott, & ideir] 
( Il975f) (CAK model). In our mo del, the line force is e valuat ed by 
using the analytical formulae of [Stevens & KallmanI ( Il990l) who 
parametrised the force in terms of the photoionization parameter, 
defined as ^ = 4n J^x/n where J^x and n are the local X-ray flux 
and the number density of the gas, respectively. This parametrisa- 
tion of the line force is computationally efficient and it provides 
good estimates for the number and opacity distribution of spectral 
lines for a given ^. 

Note that we account for the attenuation of X-ray radiation 
by computing X-ray optical depth in radial direction; however, we 
assume that the gas is optically thin in UV. The assumption of no 
attenuation of UV radiation h as to be made in our forma lism since 
the photoionization model of [Stevens & KallmanI ( ll99 Cf) does not 
include the UV attenuation; therefore, to be consistent with the their 
model. We will discuss the validity of this assumption later in Sec- 
tion|4!4l 

With the simplification above, only the corona (X-ray) radi- 
ation contributes to bringing atoms to very high ionization states. 
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We assume that the corona contributes to the radiation force due to 
electron scattering, but not to the line force although some metal 
lines in the soft X-ray may contribute to the line force in some 
cases. On the other hand, the disc radiation (UV) contributes sig- 
nificantly to the line force and the force due to electron scattering. 

In all the models presented in this work, we assume that the 
inflowing gas at the outer boundary is very slowly rotating, i.e. the 
rotation velocity of the gas is much smaller than the Keplerian ve- 
locity. We use the same specific angular momentum distribution of 
the gas at the outer boundary as in Paper III. 



2.2 Hydrodynamics 

We employ hydrodynanrical (HD) simulations of the outflow from 
accretion on to a central part of AGN in 2-D (a xisymmetric), us- 
ing the ZEUS-MP code (cf. iHaves et al.ll2006l) which is a mas- 
sive ly parallel MPI-implemented version of the ZEUS-3D code 
(cf.. lHardee & Clarkell992l : IClarkell99d) . The ZEUS-MP is a Eu- 
lerian hydrodynamics code which uses the method of finite differ- 
encing on a staggered mesh with a second -order-accurate, mono- 
tonic advection scheme jHaves et al.ll2006h . To compute the struc- 
ture and evolution of a flow irradiated by a strong continuum radi- 
ation of AGN, we solve the following set of HD equations: 
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where p, e, P and v are the mass density, energy density, pressure, 
and the velocity of gas respectively. Also, g is the gravitational 
force per unit mass. The Lagrangian/co-moving derivative is de- 
fined as D/ Dt = d/dt + v ■'V . We have introduced two new com- 
ponents to the ZEUS-MP in order to treat the gas dynamics more 
appropriate for the flow in and around AGN. The first is the acceler- 
ation due to radiative force per unit mass (g,.^^) in equation|2l and 
the second is the the effect of radiative cooling and heating, simply 
as the net cooling rate (C) in equation [3] We assume the equation 
of state to be in the form of P = (7 — 1) e where 7 is the adiabatic 
index, and 7 = 5/3 for all the models presented in this paper. Our 
numerical method used in this paper are identical, in most aspects, 
to that described in Papers I and II. Next, we briefly discuss our 
implementation of ^^.^^^j and C. 



.?^disc oc |cosS|. Consequently, the disc radiation flux at a distance 
r from the centre can be written as 



2 Icosel 



(5) 



where 6 is the polar angle (angle between the disc normal and the 
position vector r). The leading term 2 in this expression comes 
from the normalisation of the flux. 

To evaluate the r adiative acce l eration due to lin e force, we 
follow the method in IProga et alj ( I2OO0I) (see also IProga et al.l 
1998) who applied the modified version of CAK model (see also 
ICastoretal.ill975h . Their mode l works under the assumption of 
the S obolev appr oximation (e . g., Soboleviri957l : ICastoJ I97(]| ; ILucvI 
Il97lh . Following lProga et all j2000l) . the radiative acceleration of a 
unit mass at a point r can be written as 



= f [1 + M] 



a J {r,n) 



n dO, 



(6) 



where /, Q. and are the frequency-integrated continuum inten- 
sity in the direction n, the solid angle subtended by the source 
of continuum radiation, and the electron scattering cross sec- 
tion, respectively. The force multiplier is a function of optical 
dep th parameter r' which is s imilar to the Sobolev optical depth 
(c.f. lRvbicki & Hummejl978h . Th e parameters in M are func tions 
of the photoionization parameter ^ jStevens & Kallmarilll99Clll . Us- 
ing eqs. Q, ([5}, (|6]l, and the luminosity ratios (/* and /disc), the 
radiative acceleration term in equation|2]can be reduced to 



9rad 



[/* +2/disc(l + 7W) |cos9|] f. 



(7) 



To ev aluate the g a s tem perat ure, we follow the method de- 
scribed by Proga et al. (2000) and Proga & Kal lmarJ ilOOA The 
model includes some effects of photoionization. The gas is assumed 
to be optically thin to its own cooling radiation. The radiative pro- 
cesses included are Compton heating/cooling. X-ray photoioniza- 
tion and recombination, bremsstrahlung, and line cooling. The net 
cooling rate depends on the photoionization parameter ^ and the 
characteristic temperature of the X-ray radiation (Tx ). For more 
detail, readers are referred to Papers I and II (see also lProga et al.l 
(200(]|) for our implementations of g^^^ and C. 

Note that the angular distribution of the radiation in equa- 
tion|7]and also the luminosity ratios (/, and /disc) are fixed during 
the whole simulations. In other words, the radiation from the disc 
(equation [5j and the X-ray corona region (equation are treated 
as inner boundary conditions; hence, they are not self-consistently 
determined while the simulations proceed. 



2.3 Radiation force and radiative heating/cooling 

As mentioned earlier, we consider two different continuum radia- 
tion sources in our models: (1) the accretion disc, and (2) the cen- 
tral spherical corona. In the point-source approximation limit, the 
radiation flux from the central X-ray corona region can be written 

as 



(4) 



where r is the radial distance from the centre (by neglecting the 
source size). Here we neglect the geometrical obscuration of the 
corona emission by the accretion disc and vice versa. On the other 
hand, the disc radiation depends on the polar angle 8 because of 
the source geom etry. Again following Papers I and II (see also 
IProga et aljfl998h . the disc flux JTdisc is assumed to be radial and 



2.4 Coupling of Mass-Accretion Rate and Luminosity 

In our previous models (Papers I, II and III), the accretion luminos- 
ity of the system is kept fixed. Consequently, the radiative feedback 
may not be consistent with the mass supplied to the inner region of 
AGN transferred from our inner boundary. To overcome this possi- 
ble inconsistency of the accretion luminosity and the mass inflow 
rate, we relax the assumption of the constant luminosity. The lu- 
minosity is now coupled to the mass-accretion rate (A/a) in the 
following way. At each time step, the accretion luminosity is ad- 
justed based on the history of the mass accretion rate through the 
inner boundary. However, since our inner boundary is rather large 
compare to r,, we introduce a lag time (t). The lag time roughly 
corresponds to the accretion time scale for the gas to reach the 
BH from the inner boundary of our computational domain. This 
should naturally cause a self-regulation of mass-accretion rate and 
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the amount of radiative feedback, i.e. a higher mass accretion rate 
hence a higher accretion luminosity will more strongly push sur- 
rounding gas outward, and will naturally slow down the accretion 
process. 

The mass accretion rate of the BH with mass A/bh at a given 
time t is Ada, (t), the accretion luminosity is given by 

L^t^^^vGMM^ (8) 
Rs 

where rj, G and Rs are the rest mass conversion efficiency, the 
gravitational constant, and the Schwarzschild radius of the BH, re- 
spectively. The growth rate of black hole mass is negligible here in 
the time scale used in our simulation. We estimate the mass accre- 
tion rate at a given time t by computing the time average of the mass 
inflow rate (M) at the inner boundary as our simulation proceeds. 
If the time duration used for the averaging the mass accretion rate 
is At, then we can write the average mass accretion rate (Afa) as 

M.{t) = ''--^l \ . (9) 

In our models we set At = r, and the denominator simply becomes 
r. In the 'standard' a disc models, the accretion timescale (tacc) or 
the lag time r can be approximated as t^cc ~ O {tdyn/a) where 
tdyn and a are the dynamical timescal e at the outer radius (r^) o f 
the disc and the a viscosity parameter dShakura & Sunvaev|[T973h . 
In this paper, we assume = 2.6 x lO^^cm and a = 0.1, which 
give tdvn ~ lO^s, and consequently r — face ~ lO'' s (see also 
ICiotti&OstrikeJl200ll for a similar treatment). Note that the time 
scale for the free-falling gas (tff ) to reach the centre from the in- 
ner boundary of our computational domain is approximately 10* s 
which is same order of magnitude as the dynamical timescale men- 
tioned above. For the models which has a steady state solutions, we 
find that the basic flow properties (e.g. morphology, mass-accretion 
rate, and gas temperature) are rather insensitive to the value of r. 

A caveat in the method described here is that we assume all 
the inflowing gas which crosses the inner boundary reaches the ac- 
cretion disc and contributes to the accretion luminosity. In reality, 
some fraction of the gas which falls inward of the inner boundary 
may not reach the accr etion disc due to, again, strong disc radiation 
(e.g. lProga et alj2000l) . 



In the r direction, the gird is divided into 140 zones in which the 
zone size ratio is fixed at Ark+i/Ark = 1.04. 

For the initial conditions, the density and the temperature of 
gas are set uniformly, i.e. p — po and T = To everywhere in the 
computational domain (Table [T](. The initial velocity of the gas is 
assigned by assuming the same specific angular momentum distri- 
bution used in Papers II and III. At the inner and outer boundaries 
(ri and To), we apply the outflow (free-to-outflow) boundary condi- 
tions, in which the field values are extrapolated beyond the bound- 
aries using the values of t he j^host zones residing outside of normal 
computational zones (see lStone & Norman|[T99a for more details). 
At r = To, all HD quantities are assigned to the initial conditions 
(e.g. T — To and p — po) during the evolution of each model; how- 
ever, this outer boundary condition is applied only when the gas is 
inflowing at r = 7-o, i.e. when Vr < 0. The radial component of the 
velocity is allowed to float (unconstrained) when Vr > at r = Tq. 
In Paper II, we also applied these conditions to represent a steady 
flow condition at r = Tq. We found that this technique leads to 
a solution that relaxes to a steady state in both spherical and non- 
spher ical accretion with an outflow (see also IProga & Begelmanl 
l2003h . This imitates the condition in which a continuous supply of 
gas is available at r = Tq. We also consider a several cases in which 
the temperature atr — Vo is unconstrained (Section|43}. 

2.6 Reference Values 

Important reference physical quantities relevant to our systems are 
as follows. The Compton radius, Rc = GMBHpmp/kTc, is 
8 X 10^* cm or equivalently 9 x 10* r, where Tc, p and nip 
are the Compton temperature, the mean molecular weight of gas 
and the proton mass, respectively. We assume Tc = 2 x 10^ K 
and p — 1. The corresponding speed of sound at infinity is 
Coo = (ikTc/ p mpy^^ = 4 X lO^cms^^. The correspond- 
ing Bondi radius ( lBondilll95 j) is _Rb = GA/bh/c^ = 4.8 x 
10^* cm while its relation to the Compton radius is Rb ~ 7~^-Rc- 
The Bondi accretion rate (for the isothermal flow) is Mb = 
3.3 X lO^^gs"^ = O.52M0yr~\ The corresponding free-fall 
time (tfi) of gas from the Bondi radius to the inner boundary is 
2.1 X 10^^ sec = 7.0 x 10^ yr. The escape velocity from the inner 
most radius (ri = 500 r*) of the computational domain is about 

7.7 x 10* kms"\ 



2.5 Model Setup 

The following parameters are common to all the models presented 
here, and are exactly same as in Papers I, II and III. We assume 
that the central BH is non-rotating and has mass Mbh ~ 10* Mq. 
The size of the disc inner radius is assumed to be r* = 3 -Rs = 
8.8 X lO'^'^ cm. The rest mass conversion efficiency (77) in equa- 
tion [8] is assumed to be 0.0833. The Eddington luminosity of the 
Schwarzschild BH, i.e. 47rcG AfBH /o"e for our system is about 
1.3 X 10*" erg s"^ or 3.3 x 10* Lq. The fractions of the luminos- 
ity in the UV (/uv) and that in the X-ray (fx) are fixed at 0.95 and 
0.05 respectively, as in Paper I (their Run C) and in Paper II (their 
Run Cr). These two parameters determine the shape of the under- 
lying continuum emission of the central source (cf. Section lZ2t . 

The following ranges of the coordinates are adopted: ri ^ r ^ 
ro, < 6* ^ 7r/2 where n = 500 r* and To = 2.5 x 10^ r,. The 
radius of the central and spherical X-ray corona region r, coincides 
the inner radius of the the accretion disc. In our simulations, the 
polar angle range is divided into 64 zones, and are equally spaced. 



3 RESULTS 

We examine the effect of changing the outer boundary conditions 
on the gas property, flow morphology and gas dynamics. Here, we 
focus on two key parameters, i.e. density (po) and temperature (To) 
at r = To. We consider three different temperatures: To = 2 x 
10*^, 2 X lO'^ and 2 x 10* K with a wide range of p^. In total, 22 
combinations of po and To are considered. The model parameters 
along basic results are summarised in Table. [T] 

3.1 Dependency of Mass Outflow Rates on Accretion 
Luminosity 

First, we examine how mass outflow rates of systems depend on 
the accretion luminosity that is self-consistently determined. Fig.[T] 
shows the mass outflow rate at the outer boundary Afout (''o) plot- 
ted as a function of the accretion luminosity, or equivalently as 
a function of the Eddington ratio F (= L/LecW where LEdd = 
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Table 1. Summary of models with different combinations of po and To 





po 


To 


Min (n) 


r 


A'/out (r-o) 


Variability 


Outflow Morphology 


Model 


(10 gem 


(10' K) 


(1025 gs ^) 




(10^^ gs '-) 






1 


1 


0.2 


5 1 (n 002)* 


0.30(0. OOl)'" 


OfO 01* 


steady 


no outflow 


2 


2 


0.2 


8 6fO 141 


51 fO 01 1 


1 4fl 31 


steady 




3 


3 


0.2 


11(0.36) 


0.65(0.02) 


3.4(2.1) 


semi-steady 


wide polar 


4 


5 


0.2 


15(1.8) 


0.89(0.11) 


6.2(3.1) 


semi- steady 


wide polar 


5 


10 


0.2 


22(2.3) 


1.3(0.1) 


13(4.6) 


semi-steady 


wide polar 




20 


0.2 


36('3 8") 


2.1(0.2) 


71 


semi-steady 


disc wind 


7 


50 


0.2 


63(8.6) 


3.8(0.5) 


21(16) 


semi-steady 


disc wind 


g 


100 


0.2 


93(6.2) 


5.5(0.4) 


18(13) 


semi-steady 


disc wind 


9 


200 


0.2 


130(5 01 


7.7(0.3) 


20(12) 


semi- steady 


disc wind 


10 


0.4 


2 


2.4(0.01) 


\A(0 01 1 


OfO 01 


steady 


no outflow 


11 


1 


2 


8.3(0.36) 


0.49(0.02) 


2.0(2.3) 


semi- steady 


narrow conic 


12 


2 


2 


15(6.5) 


0.89(0.38) 


9.4(11) 


semi-steady 


conic 


13 


4 


2 


26(13) 


1.5(0.8) 


24(22) 


semi-steady 


wide conic 


14 


10 


2 


56(31) 


3 3('l 81 


66(91) 


semi- steady 


wide conic 


15 


0.5 


20 


1.3(4.9) 


0.074(0.292) 


3.4(4.5) 


spiky 


almost spherical 


16 


1 


20 


2.8(5.6) 


0.17(0.33) 


4.3(4.2) 


spiky 


spherical/polar 


17 


1.5 


20 


7.5(13) 


0.45(0.77) 


7.1(9.1) 


spiky 


polar 


18 


2 


20 


19(0.072) 


1.1(0.043) 


13(0.51) 


steady 


polar 


19 


3 


20 


29(0.080) 


1.7(0.047) 


31(2.6) 


steady 


polar 


20 


4 


20 


39(0.31) 


2.3(0.02) 


54(16) 


steady 


conic 


21 


6 


20 


58(3.0) 


3.5(0.2) 


113(47) 


semi-steady 


conic 


22 


8 


20 


77(3.6) 


4.6(0.2) 


169(66) 


semi-steady 


wide conic 


23* 


1 


2 


8.3(0.38) 


0.49(0.02) 


2.0(1.8) 


semi-steady 


nan'ow conic 


24* 


2 


2 


16(4.7) 


0.95(0.28) 


12(7.0) 


semi-steady 


conic 


25* 


4 


2 


40(13) 


2.4**(0.77) 


14(11) 


semi- steady 


conic 


26* 


10 


2 


140(91) 


8.3**(5.4) 


23(41) 


semi- steady 


conic 


27t 


1 


2 


4.7(0.29) 


0.27(0.02) 


5.4(1.5) 


semi-steady 


conic 



(*) Luminosity is Umited to F sC 1 and not completely self-consistent. 

(**) r value corresponds to the luminosity which would have been achieved in the system is if the luminosity is computed from the mass accretion rate listed 

in the previous column. 

(t) Fixed luminosity model in lProga etalj )2008h (their Run Cr) and lKurosawa & Progal )2009l) (their Model II). 
(I) Values in brackets are the standard deviations of the time series values. 



Table 2. Summary of power-law fits for the Mout- F relations 



Main wind driving force T F range g** 

(10^ K) 



Radiation 


all* 


~ 0.5 to - 


- 5 


2.0 (±0.1) 


Radiation 


0.2 


~ 0.5 to ~ 


1.4 


2.3 (±0.3) 


Radiation 


2 


~ 0.5 to - 


- 3 


1.8 (±0.2) 


Radiation 


20 


~ 0.9 to - 


- 5 


1.8 (±0.1) 


Thermal 


20 


< 0.7 




0.4 (±0.1) 


Radiation ('disc-wind') 


0.2 


>1.4 




0.1 (±0.1) 



(*) Used the points from all three temperatures (excluding points for Mod- 
els 6, 7, 8, 9, 15, 16 and 17) 

(**) The index of th e power-law, i.e. M put oc F'. The radiation-driven 
stellar wind model of lCastor et al.l jl975h finds q = 1.67. 



AncGM^Yi/cre ) for different combinations of po and To. The ac- 
tual values for Mout (^o) and F (along with the corresponding mass 
inflow fluxes at the inner boundary A/in) are also listed in Table. [T] 
Note that Afout values used here are the time-averaged values since 
the mass inflow rate A/in or equivalently mass accretion rate A/a in 



some of the models (as also indicated in the table) show some de- 
gree of variability, which will be discussed later in Section [T2l 

The figure shows that most of the data points follow the lines 
that cross the panel almost diagonally, displaying strong correla- 
tions between A/out and F. The outflows formed in the models with 
r > 0.5 are mainly mdiatively driven whereas those in the highest 
temperature cases (To = 2 x 10* K) with F < 0.5 (Models 15, 16 
and 17) are mainly thermally driven. 

The models with the lowest temperature (To = 2 x 10® K) 
follow the main trend defined by the data points for the models 
with the radiation-driving dominated outflows for 0.5 < F < 1.5. 
No outflow forms for F < 0.5. The mass outflow rates become 
saturated at higher F values (F > 1.5), i.e. the strength of the 
correlation weakens dramatically. The cause of the saturation seen 
in the the low temperature models will be discussed later in Sec- 
tion 13.31 The models with To = 2 x 10^ K shows a correlation 
between A/out and F for all models with F > 0.5. Again no out- 
flow forms for F < 0.5. The figure also shows a point for the 
model run (Model 27) with the same parameters as in Model 1 1 
(po = 1 X 10"2i gcm"^ and To = 2 x lO'^ K), but with a fixed 
accretion luminosity (F = 0.6). Note that this model is equivalent 
to Run Cr in Paper II, and Model II in Paper III. These fixed F (con- 
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Figure 1. Comparisons of the relations between the mass outflow rates 
(A/out through the outer boundary) and the Eddington ratio (F) for the 
models with three different ambient gas temperatures (To) : 2 X 10^ (cir- 
cles), 2 X 10^ (squires) and 2 X 10* K (diamonds). The parameters used 
for the models/data points are summarised in Table [T] The points along 
the diagonal lines (thick lines) clearly shows that Mout correlates with F 
for all To. The strong correlation is seen in the lowest temperature mod- 
els (To = 2 X 10«K) for 0.5 < F < 1.5, but J\/out becomes satu- 
rated at the higher F values (F > 1.5; cf. Sec. l3.3t . For the models with 
To = 2 X 10^ K, a correlation is seen for all the points with F > 0.5. The 
high temperature models (To = 2 X 10® K) show two distinctive regions 
in which the points follow two different slopes. The outflows in the models 
for F < 1 (Models 15, 16 and 17) are mainly thermally driven, but the 
models with F > 1 are mainly radiation-driven. The thermal winds found 
here are artifacts of the itmer boundary. A simple power law is used to fit 
the points for the models with the outflows mainly driven by radiation, for 
each temperature: 2 X 10^ (thick dashed line), 2 X 10^ (thick solid line) 
and 2 X 10** K (thick dash-dot line). A separate power law (thin dashed 
line) is used for the points with four largest F (Models 6, 7, 8 and 9) in the 
To = 2 X 10'' K models, and another power law is used for the models 
with the thermally driven outflows with To = 2 X 10* K (Models 15, 16 
and 17). The indices of the power law fits are summarised in Table|2] The 
mass-lo ss rate predicted by the radiation-driven stellar wind model (dotted 
line) bv lCastor et al.l i 19751) (CAK model) is also shown for a comparison. 
The arrows near the bottom of the panel indicate the approximate location of 
the Eddington ratios for the models below which an outflow does not form. 
The solid and dashed arrows are for the models with To = 2 X 10^ and 
2 X 10^ K, respectively. The constant luminosity (F = 0.6) model (star), as 



in Papers II and III (with To = 2 X lO'^ K and po 
is also shown for a comparison (Model 27). 



1 X 10^21 



gem ^), 



slant L) models are located near the main trend of the A/out - F re- 
lation found here. This reassures that our previous models parame- 
ters (F, Po and To) were reasonable, and were somewhat consistent 
with each others. The high temperature models (To — 2 x 10** K) 
show two distinctive regions in which the A^out - F relation follows 
two different slopes. As mentioned earlier, the outflows found in the 
models with F < 1 (Models 15, 16 and 17) are mainly thermally 
driven, but the models with F > 1 are mainly radiation-driven. The 
change in the main wind driving mechanism (from thermal pres- 
sure to radiation pressure) is the cause of the change in the strength 
or the slope of the A/out - F coiTelation. The cause of the super- 



Eddington accretion luminosity (F > 1) will be discussed later in 
Sectionlim 

A simple power law (A/out oc F') is use to fit the points for the 
models with radiation-driving dominated outflows, separately for 
each temperature. We find the power law indices (q) as 2.3 (±0.3), 
1.8 (±0.2) and 1.8 (±0.1) for To = 2 x 10^ 2 x lO'' and 2 x 10* K, 
respectively. These values are very similar to each others, i.e. the 
slope of the correlation is independent of or very insensitive to the 
outer boundary temperature. The slope for To = 2 x 10® K seems 
to be slightly higher than those of the higher temperature models; 
however, when the point corresponding to Model 2 (po = 2 x 
10^^^ gem"'') is excluded from the line fit, we obtain the power 
law index g = 1.9 which is in agreement with those of the models 
with the higher temperatures. We note that the outflow in Model 2 is 
relatively weak and it may not be entirely radiation-driven as its F 
is very close to that of Model 1, which does not produce an outflow. 

The figure also shows th e mass-loss r ate pr edicted by the 
radiation-driven wind model of ICastor et al.l ( 1 19751) (CAK model) 
for a comparison. The CAK model predicts q = 1.67 for F ^ 1 
and g J? 2 for F ^ 1 which are very similar to the slopes found 
here. Their model assumes spherical symmetry, and the mass out- 
flow rates becomes infinity as F approaches to 1. The differences 
between their model and our model are not only due to the assumed 
symmetries (spherical versus axial), but also in the implementa- 
tion of the underlining luminosities. Our wind driving luminosity 
is coupled to the mass inflow/accretion rate of the system while 
their luminosity is fixed constant and they do not have inflows in 
their model. 

Motivated by the insensitivity of the slope of A'/out - F curves 
on To, we compute the slope using the points from all three tem- 
peratures (excluding the thermally-driving dominated models and 
the lowest temperature models with F > 1.5: Models 6, 7, 8, 9, 
15, 16 and 17). The corresponding slope is found as 2.0 (±0.1). 
A separate power law is applied to the points with four highest F 
(Models 6, 7, 8 and 9) in the To = 2 x lO" K models, and its 
index is 0.1 (±0.1). For the models with the thermally driven out- 
flows with To = 2 x 10*K (Models 15, 16 and 17), the power 
law index is 0.4 (±0.1). Tabled summarises the power law indices 
found in this analysis. We note that theoretical a nd numerical stud- 
ies of radiation driven stellar and disc win ds (e.g. iCastor et aljl975l ; 
IProga, Stone. & Drew|[l993 ; |ProgJ 19991) predict q^ 1.67, which 
is quite close to the values found in this study. 

When the outer boundary density po becomes too small (po < 
1 X 10"^^ g cm""'), as in Model 1, no outflow is formed for 
To = 2 X 10® K. Similarly, no outflow is found for the models with 
Po < 4 X 10"^^ gcm-^* as in Model 10, for To = 2 x lO'^ K. 
Unlike in the lower temperature models, an outflow forms for 
To = 2x10* K, even with a relatively small outer boundary density 
(po ~ 0.5 x 10~^^ gem""'). The minimum outer boundary den- 
sity Po"" or the minimum Eddington ratio p""" above which an 
outflow can be formed depends on the ambient temperature (To). 
The higher the value of To, the smaller the value of r™'". For ex- 
ample, F™'" ~ 0.3, ~ 0.2 and < 0.05 for To = 2 x 10®, 2 x lO'^ 
and 2 x 10* K, respectively. When the outer boundary temperature 
becomes higher, the average temperature of the gas in the tempera- 
ture becomes also higher. This makes the gas to be driven thermally 
more easily; hence, the outflow starts to form even at lower value 
of F. Note that the Eddington ratio F, by definition, only considers 
the radiation force due to Thomson scattering. The radiation force 
considered here includes also line scattering process (see e.g. Pa- 
per I for our implementation of the radiation force); hence, the total 
radiation force could exceed the gravity even for F < 1. 



The correlation between Afout and F (for the radiation driven 
outflow cases) may continue for higher luminosities (i.e. F > 5), 
but we are not be able to confirm this due to a limitation of our 
code. Strong shocks form in the flows in the models with F > 5 
(for To = 2 X 10'' and 2 x 10* K models), and the code fails to 
handle them. Unfortunately, this prevents us from exploring models 
with very high accretion luminosities. 

3.2 Time-Dependent Behaviours of Mass-Accretion Rates 

Not all the models we consider here reach steady state solutions 
and have smooth inflows (and outflows). Some models do show 
variability in the mass-accretion rate (Afa) (the mass flux across the 
inner boundary of our computational domain). To quantify the de- 
grees of variability, we compute the standard deviations of the time 
series values of and F for each model. The results are shown 
in Table [T] Further, to compare the amount of variability from dif- 
ferent models more directly, we normalise the standard deviation 
by the time-averaged values of Afa. In other words, we compute 
the ratio ((t„) of the standard deviation of Afa to the time-averaged 
values of Afa. The distribution of the normalised standard deviation 
an is shown in Fig.|2] The figure shows an as a function of F for 
Models 1-22 in TableQ 

We categorise the time-dependent natures of in our sim- 
ulations into three types: (1) steady, (2) semi-steady and (3) spiky, 
based on the value an- We define a model as steady when an < 
0.01, and as semi-steady when 0.01 < (t„ < 1.0. Further, it is 
defined as spiky when cr„ > 1.0. Based on this definition, the vari- 
ability type of each model is assigned and placed also in Table [T] 
Out of 22 models considered here, 6 are found as steady, and 13 
as semi-steady. Only 3 models are found as spiky. Fig.[3]shows ex- 
amples of the three different types of A/a variability found in our 
simulations. 

In the models with steady A/a, inflows and outflows are 
smooth and the level of A/a variability (the change with respect 
to the time averaged value of A/a) is typically less than 1 per cent. 
In the models with semi-steady A/a, small-sized and relatively high 
density structures are occasionally formed in the flows, and some 
of them fall toward the centre. This causes moderate levels of vari- 
ability (<200 per cent) in a relatively short timescale (~ 10^^ s). 
The models marked as spiky occurs only in the high temperature 
models (To = 2 x 10** K), and in the cases when the outflows is 
thermally driven. The relatively high temperature of gas keeps the 
sound speed to be high. As a result, the outflows in these mod- 
els are sub-sonic. The effect of a small density perturbation in the 
flow that occurs near the inner boundary propagates though the en- 
tire computational space in relatively small timescale (because of 
the high sound speed) compared to the dynamical timescale of the 
flows. The acoustic wave generated near inner boundary propagates 
outward and bounces back at r = ro in the hot media, when the 
outflow is relatively weak and when the radiative force is relatively 
weak. This causes rather stochastic sharp peaks and dips (spikes) 
in the A/a variability curves. The level of variability is rather large 
(> 100 per cent) in these models. We do not consider the models 
with spiky A/a (i.e. Models 15, 16 and 17) as physical because they 
are affected by the inner boundary. The inner boundary is too far 
from the BH for the inflow to become supersonic in those cases; 
hence, it can reflect the flows and turn inflows into outflows. 

When we exclude the unphysical spiky A/a models (i.e. Mod- 
els 15, 16 and 17) from Fig.lS we find a very weak correlation be- 
tween an and F, i.e. the larger the accretion luminosity, the amount 
of variability tends to be larger. The distribution of an seems to be 
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Figure 2. Distribution of the normalised standard deviation (t„ of the mass 
accretion rates computed at the inner boundary for Models 1-22 in Table[T] 
is shown as a function of the Eddington ration P. models are categorised 
as steady, semi-steady or spiky when (t„ < 0.01, 0.01 < cr„ < 1.0 
and (T„ > 1.0, respectively. The boundaries between different types ai'e 
indicated by horizontal lines {dashed line). The models which belong to the 
.steady and semi-steady types are indicated by circles, and those belong to 
the spiky are indicated by squares. A very weak correlation between (j„ and 
r is seen for the data points belong to the steady and semi-steady types. 

continuous (as a function F). There is no clear separation of the 
two populations: the steady and semi-steady types, according to the 
figure. 

3.3 Disc Wind Like Solutions 

In most of the models explored (Table [T), the flow morphology 
are quite similar to the ones presented earlier (Papers I, II and III) 
e.g. bipolar outflows and fragmented dense cold cloud outflows. 
Here we report a new morphological type of solutions which has 
not been seen in our previous models. Among the low temperature 
(To = 2 X 10^) models, when the density at r = To is relatively 
high (po > 5.0 X 10"'^'^ g cm~^), the flow starts to resemble a 
radiation driven disc wind which arises from the very inn er part 
of accretion discs (e.g. lProga et al.|[l998l : |Proga et al.lll999l) , and a 
ther mally-driven disc wind which arises from re latively large radii 
(e.g. lWoods et al .111 9961 : [Proga & Kallmarill2002h . 

Fig Ejdemonstrates how the wind configuration transforms to 
a disc wind like solution. For this purpose, we compare the gas 
flows from Models 5, 8 and 9 which have the same outer bound- 
ary temperature To = 2 x 10® K but have different values of po, 
i.e. 1 X 10"^", 1 X 10"^^ and 2 x lO"^** gcm"^ , respectively. 
In Model 5, which has the lowest po value among the three, the 
high density gas streams falling toward the centre in relatively wide 
range of polar angles (50° ^ 6 < 90°). The outflow is formed 
only in the mid polar angles 30° ^ < 50° . When po is increased 
10 times (Model 8), the infalling dense streams of gas seen in the 
previous model are converged toward equatorial plane, and form 
a relatively dense disc-like structure. The outflow now occurs at 
larger polar angles and in wider angle range (45° < S < 80°). 
The outflow resembles a disc wind. However, the high density gas 
near r = ro at lower polar angles 6 < 45° prevents the gas from 
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Figure 3. Examples of different types of time-dependent mass-accretion 
rate variability patterns. The panels show the mass-accretion rates as a func- 
tion of time for selected samples of models which show steady (top), semi- 
steady (middle) and spiky (bottom) variabilities (see Table[T). The number 
on the upper right hand comer of each panel indicates the con'esponding 
model number in Table[T] The mass accretion rates on the vertical axes ai'e 
in units of the time-averaged mass accretion rate of each model which are 
39, 15 and 7.5x10^^ gs^^ from the top to bottom panels, respectively. 



forming outflows and keeps the gas relatively hot in the polar direc- 
tion. In this case, no fast low- density polar wind is formed unlike 
the ones seen in the models of IProga et ai] ( Il998l) . This difference 
is cau sed by the difference in the treatment of boundary conditions, 
i.e. in lProga eTal, ( 1998). the density and temperature at r = ro are 
unconstrained. When we increase the density po by a factor of 2 
(Model 9), the radiative force becomes strong enough to push away 
the high density gas in the polar direction, and the fast low-density 
polar wind is formed. However, in this case, the outflowing mass- 
flux in the polar direction is negligible (due to its very low den- 
sity) compared to the outflow in the higher polar angles {6 > 45°). 
The accretion region is now converged to a much thinner equatorial 
disc, and a very wide outflow (45° < ^ < 85°) or a disc wind is 
formed. 

The saturation of the mass outflow rates at the high F val- 
ues for To = 2 x 10® K models seen in Fig. [T] (Section [Z4) can 
be explained qualitatively by the change of the outflow pattern to a 
disc wind like configuration. Before the saturation occurs, the polar 
outflow can become wider by squeezing the high density inflowing 
stream toward the equatorial regions. This allows the mass outflow 
rate to grow proportionally with the mass-accretion rate or with the 
accretion luminosity F. In the higher F models, e.g. Models 8 and 
9, the dense inflowing gas streams are converged all the way to the 
equatorial plane (9 ~ 90°); hence, the polar angle range, in which 
outflow occurs, can not be increased any more. This limit in the size 
of the outflow polar angle range is the cause of the saturation of the 
mass outflow rates (the flat part of the A/a-F relation in Fig. [TJ. 



In addition, once the flow becomes disc-wind like, it is harder for 
a system to produce stronger outflows for the following reason. In 
our models, the radiation peaks in polar directions (cf. cos 6 depen- 
dency of the radiation flux as in Sections |2. 1 1 and \4~]\ : hence, the 
most of the radiation escapes in the polar direction. In the disc-wind 
like flow morphology, however, there is almost no gas to be blown 
away by the radiation in the polar regions (Fig. |4ll. This results in 
the insensitivity of the mass outflow rate with F. Near the equatorial 
plane where most of the infalling gas is concentrated, the radiation 
is very weak because of the cos 9 dependency of the disc radiation, 
and because of the very high optical depth of the disc like structure. 
Finally, the temperature map for the disc-wind like solution shows 
that the temperature of the wind is relatively high (T ~ 10^ K). 
This temperature is too high for the radiation force due to line pro- 
cesses (line force hereafter) to be effective. We confirmed that the 
main wind driving force for the wind in this case is the continuum 
radiation force (due to electron scattering). 



4 DISCUSSIONS 

4.1 Accretion with Super-Eddington Luminosity 

Many of the self-consistently determined accretion luminosities in 
the models presented in Section \3A\ (see Table [T] and Fig.O are 
found to be above the Eddington luminosity (F > 1). At first, this 
may be puzzling since one would expect that the high luminosity 
hence strong radiation pressure induced by the high mass accretion 
rate should regulate the amount of luminosity and keep F < 1. 
This might be the case if the gas were not rotating and the radi- 
ation were spherically symmetric; however, the situation changes 
when the source geometry of the strong UV radiation is not spher- 
ical. Because of the disc geometry, the flux of the disc radiation 
depends on the polar angle. It radiates more in the polar directions 
(J?-dis(. oc |cos^j), as mentioned in Section \22\ and it turns in- 
flows into outflows. To see how this assumption leads to the super- 
Eddington accretion solutions found in our models, consider the 
ratios (Fg) of the radial forces per unit mass (accelerations) by ra- 
diation (flrad) and gravity (g) (cf. equation|2]l. To simplify our argu- 
ment, we consider only the radiation force due to electron scatter- 
ing and exclude that due to spectral line processes here. By setting 
= in equationl?] the radiative acceleration can be written as 



.9rad — 



L 



if* + 2/dis 



(10) 



where and L are the electron scattering cross section and the 
total accretion luminosity respectively. The superscript e on the left 
term indicates that the acceleration is only due to electron scat- 
tering. The fractional luminosities of the X-ray emitting spherical 
region and that of the UV emitting accretion are assumed to be 
/, = 0.05 and /disc ~ 0.95 (Section [2.1t . Using equation [Tol and 
the gravitational acceleration g = —G Msn/r^ (assuming it is 
Newtonian and a point-mass), the ratio becomes 



ffrad 



r(/, +2/disc cos^) 



(11) 



where F = (7eL/(47r (?A/bh) and is the Eddington ratio. Fig. |5] 
shows the values of Fg plotted as a function of the polar angle 9 for 
a range of the Eddington ratios (0.1 < F < 10). We consider only 
the radiation force due to electron scattering. Therefore, the lines 
shown in the figure are only lower limits of Fg values for at a given 
9 angle. The figure shows Tg < 1 (the outward radiation force is 
smaller than the inward gravitation force) at all 9 for F < 0.5 cases. 
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Figure 4. A transition to a disc-wind like solution in the low temperature models (To = 2 X 10® K). The density (in logaiithmic scale) over-plotted with the 
directions of poloidal velocity as arrows (left column) for Models 5 (fop). 8 (middle) and 9 (bottom) are shown. The figures are placed in order of increasing 
density at r = ro (po), from the top to bottom (cf. Table[T). The corresponding temperature maps (in logarithmic scale) of each model are placed in the right 
column. The colour tables for the density and temperature maps are labelled with the logarithmic values in gcm""^ and K, respectively. The top density plot 
shows the high density gas streams falling toward the centre in the wide range of polar angles (8 > 50° ), and the outflow is formed only in the mid polar 
angles 30° ^6 ^ 50°. The middle density plot shows the infalling dense gas forming a relative ly thick disc-like structure, and the outflow occurs in larger 
and wider polar angles (45° ^ 6 ^ 80°) which resembles a disc wind (e.g. IProga et aljri998l) . The high density gas near r = ro at lower polar angles 
(8 < 45°) prevents the gas from forming outflows in the polar direction. The bottom density plot shows the accretion region is now squeezed to a much 
thinner equatorial disc, and a very wide outflows (45° ^8 < 85°) or a disc wind is formed. 
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For the models with F > 1, the figure shows that Te can still be less 
than 1 at larger values of 9. For example, for F = 5 model, Ve is 
still smaller than 1 for 61 > 85° , and for F = 10 case, f > 88° . This 
critical angle Scrit at which the radiation force equals the gravity 
(Fe = 1) becomes lager as the accretion luminosity or F increases. 
The inflow of gas is potentially possible only in the polar angle 
range which has Te < 1. Therefore, the location of the accretion 
streams becomes closer and closer to the equatorial plane. The most 
important point of Fig.[5]for us is that Te can be still less than I near 
the equatorial plane which allows an accretion of gas in a super- 
Eddington luminosity environment. This is a direct consequence of 
the disc geometry which allows most of the radiation to escape in 
polar direction, but not in the equatorial direction. In addition, the 
radiative cooling in the dense gas confined near the equatorial plane 
is very efficient, and the gas temperature is kept low (e.g. see the 
bottom panels of Fig.|4l(. This also helps to accrete the matter in the 
equatorial plane. 

Interestingly, a recent 3-D radiation-hydrodynamic simula- 
tio n of the format i on of a massive (stellar) binary by accretion 
bv lKrumholz et al.l ilOO^ demonstrated that the radiation pressure 
would not be able to halt accretion totally even at F > 1 because 
of gravitational and Rayleigh Taylor instabilities which channel the 
accretion of matter on to the binary from the accretion disc, and 
the self-shielding of radiation by high density gas filaments which 
could still experience a net inward force. Note that this type of 
self- shielding does not occur in our model since the UV radiation, 
which causes 95 per cent or more of the total radiation force in our 
model, is assumed to be optically thin in our flows. 

4.2 Effect of Luminosity Limit 

The self-consistent accretion luminosity models presented in Sec- 
tion 13.11 do not limit their luminosities to be below the Ed- 
dington luminosity (i.e. F < 1) as explained above. A super- 
Eddingto n or super-critical accretion on the SMBH may be pos- 
. Jaroszvnski et aljri98(ll:lAbramowicz. Calvani, & Nobili 



sible (e.j 

ll98d : lBegelmarJ|2002l:IOhsuga et alfcooaiOhsuga et aL2005. : see 



also a review bv lAbramowica 20051 .). However, we can not be cer- 
tain whether the accretion disc would be still stable and would ac- 
crete at a such high Eddington ratio (e.g. F — 3.3 for Model 14) 
because we do not model the dynamical evolution of the accretion 
disc itself (cf. Section [231 l. When the accretion disc becomes op- 
tically thick, the radiation field of the interior of the disc itself is 
expected to become very isotropic. In general, a super-Eddington 
accretion in such environment is much more difficult than the one 
we considered in this paper. When the radiation field is isotropic, 
one would expect the accreting gas to be turned around when F > 1 
and the disc becomes optically thick; hence, it would tend to keep 
the mass-accreting rate below the super-Eddington rate. For this 
reason, we examine the possibility that the accretion luminosity is 
limited to the Eddington luminosity, and study the consequences of 
such a limit. 

We recomputed the models with the self-consistent evaluation 
of the accretion luminosity, but now the upper limit of the accre- 
tion luminosity is set to the Eddington luminosity. Fig.|6]shows the 
mass outflow rates at r = r-o plotted against the mass accretion 
rates at the inner boundary for the models with To = 2 x 10^ K 
with and without the upper limit on the luminosity. Only the models 
with high mass accretion rates (near and above the mass accretion 
rate corresponding to F = 1) are shown in the figure (cf. Table[T]l. 
The figure is similar to Fig. [T] but the mass accretion rate is used 
for the horizontal axis instead of F. The figure shows that mass out- 
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Figure 5. Ratios (Fg) of the radiation force due to electron scattering (g'^^^) 
to the gravitational force (g) as a function of spatial polar angle 8 for dif- 
ferent values of the Eddington ratios F (solid lines). The lines indicate the 
lower limit of the ratio because the radiation force due to line process are 
not included here. The horizontal dashed line indicates Fg = 1 at which 
the magnitudes of radiation force and the gravitational force are equal. In 
the polai' angle ranges below the Fg = 1 line, an inflow of gas is possible. 
This range decreases as F becomes larger, indicating that an inflow can oc- 
cur in larger polar angles, i.e. closer to the equatorial planes. Note that even 
for r = 5, Fg is still smaller than 1 for S > 85°, and for F = 10 case, 
9 > 88°. The fractional luminosities of the spherical X-ray emission and 
that of the disc UV emission are assumed to be 0.05 and 0.95 respectively 
(cf. Section im . 



flow weakens (the mass outflow rate decreases) for the models with 
larger mass accretion rates (larger than the value corresponding to 
F = 1). On the other hand, the models with smaller mass accre- 
tion rates (below F = 1 line in the figure) have nearly identical 
mass outflow rates. This behaviour is expected from the luminosity 
limiting model. The models with higher mass accretion rates are 
affected by the luminosity limit hence have smaller luminosities 
compared to the cases without the limit. A smaller luminosity nat- 
urally leads to a lower mass outflow rate. Interestingly this decrease 
or the saturation of the outflow strength (the mass outflow rates) is 
very similar to that one seen in the high F end of the low tempera- 
ture models (To = 2 x lO" K) in Section[3T](see Fig.© although 
the physical reason for the cause of the saturation is different. 



4.3 Effect of unconstraining the outer boundary temperature 

Motivated by the insensitivity of the slopes of A/out-F relations 
on the outer boundary temperature (To) found in Section [3T| we 
now investigate the effect of allowing the temperature at ro to 
be computed self-consistently instead of being fixed at some pre- 
determined value. The models presented so far have assumed that 
the gas which surrounds our computational domain is 'Comp- 
tonised' at a constant temperature To- Here, we consider 7 addi- 
tional models in which we do not constrain the temperatures at 
r = To- We now compute the temperatures self-consistently by 
solving the equation of the energy (i.e., equation O at To, as it is 
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Figure 6. Effect of limiting accretion luminosity on the mass outflow rates 
for the models with the ambient temperature To = 2 X 10^ K. The mass 
outflow rates are plotted against the mass inflow rates for the models with 
(triangles) and without (squares) the threshold accretion luminosity for a 
comparison. For the models with the luminosity limit, the maximum al- 
lowed luminosity is set to the Eddington luminosity, i.e. F = 1. The mass 
inflow rate that con'esponds to F = 1 is indicated by a vertical dotted line. 



done for any other points in our computational domain. The basic 
model parameters and results are summarised in Table[3] 

Fig. |7] shows the mass outflow rates (Mout through the outer 
boundary) plotted as a function of the Eddington ration (F) for the 
new runs, along with those of the models with To = 2 x 10'' K 
from Section [3n (see also Table [T). Both sets of the models show 
a very similar dependency of Mout on F. This is mainly because 
the average temperatures (self-consistently determined) at ro for 
the new models considered here are relatively low and similar to 
that of the fixed outer boundary temperature models i.e. To = 2 x 
10^ K. Table [3] lists the average outer boundary temperatures for 
the individual runs. 

A strong correlation between Mout and F is found for the 
models with F < 2. The new set of models also show a saturation 
of A/out beyond F > 2 as in the To = 2 x lO'^ K (fixed) models, 
and exhibit the disc wind like solutions (cf. Section |33] |. We ap- 
plied a simple power-law fit for the points with F < 2 (Models 28, 
29, 30 and 31 in Table [U, and found its index as q = 2.2 ± 0.4 
which is very similar to that of the To = 2 x 10* K models and that 
of the CAK stellar wind model (cf. Table|2j. No outflow was found 
for the models with F < 0.3. 

The self-consistently determined temperatures at r = ro in 
the new models depend on the polar angle. In general, the temper- 
ature decreases from ~ 10^ K in the polar direction (0 — 0) to 
~ lO'' K as increases toward the equatorial plane. At the higher 
polar angles (6 > 50°), relatively high density accretion streams 
are present. The high density columns of accreting gas have high 
optical depths; hence, they are very efficient for shielding the strong 
ionization radiation from the central source. This causes the lower 
temperatures at r = ro for the higher 6 values. 

In summary, the inflow and outflow properties for the new 
models are very similar to those of the low outer boundary tem- 



r, 10 



10 




A 


T unconstrained 


o 


T__=2x10'k 




CAK model 



0.1 



10 



Figure 7. Comparison of models with (circles) and without (triangles) fix- 
ing the temperature at the outer boundary (To). The mass outflow rates 
(Mout through the outer boundary) are plotted as a function of the Ed- 
dington ration (F) for the models with To = 2 X 10^ K from Section IbD 
(see Table[T) and for the new models without To constrained (see Table[3). 
Both sets of the models show a very similar dependency of Afout on F. The 
models shows strong coiTelations between Mout and F when F < 2. The 
new set of models also show a saturation of A/out beyond F > 2 as in the 
To = 2 X 10*' K (fixed) models, and exhibit the disc wind like solutions 
(see Section [33). A simple power law is used to fit the points with F < 2 
(Models 28, 29, 30 and 31 in Table|3) which show a strong correlation be- 
tween A/out and F. The index of the power law fit for the new models is 
g = 2.2 ± 0.4 which is very similar to fliat of the To = 2 X 10"^ K models 
(see Table 12). The mass-loss rate predicted by the radiation-driven stellar 
wind model (CAK model) is also shown for a comparison (dash-dot line). 
The vertical line at F 0.3 (dotted line) indicates the approximate F value 
below which a model does not produce outflows. 



perature models (To = 2 x 10^ K). The disc wind like solutions 
also appear for the models with high accretion rates (F > 2). 



4.4 Assumption of flows with optically-thin UV radiation 



As briefly mentioned in Section IXTl our model takes into account 
for the attenuation of X-ray radiation by computing X-ray optical 
depth when computing the radiative cooling/heating rate and radi- 
ation forces; however, we assume that the flow is optically thin in 
UV. The attenuation of the UV radiation is not considered here be- 
cause our model uses the method described bv lStevens & Kallmanl 
( Il99(lh who parametrised the line force by the photoionization pa- 
rameter ^. This simplification speeds up our calculation tremen- 
dously since ^ can be estimated rather quickly without comput- 
ing the detail ionization levels of the gas. The formulation of 
IStevens & Kallmanl ( IT990I) assumes that the UV radiation is opti- 
cally thin; therefore, to be consistent with their original model, our 
model must make the same assumption. 

We have checked the radial UV optical depths (between the in- 
ner and outer boundaries) of a several representative models (espe- 
cially for those with high mass accretion rates, F > 1) at a few time 
slices. For each case, the optical depths (r) have been computed 
along each polar angle (9). We found that r in general increases as 6 
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Table 3. Summary of models without To constrained 



Model 


Po 


Mn in) 

(10^^ gs-i) 


r 


A'/out (^o) 

(10^5 gs-1) 


T* 

avG 

(10^ K) 


Variability 


Outflow Morpholo 


28 


10 


8.9(0.68)* 


0.52(0.04)* 


1.0(0.032)* 


1.4 


semi-steady 


wide polar 


29 


20 


12(0.58) 


0.71(0.034) 


3.5(0.52) 


1.4 


semi-steady 


wide polar 


30 


40 


18(1.1) 


1.1(0.1) 


7.1(1.4) 


3.6 


semi-steady 


wide polar 


31 


80 


25(2.6) 


1.4(0.2) 


10(2.3) 


8.0 


semi-steady 


wide polar 


32 


160 


36(4.9) 


2.1(0.3) 


9.9(2.6) 


9.8 


semi-steady 


disc-wind like 


33 


320 


52(2.6) 


3.1(0.2) 


11(0.48) 


8.5 


semi-steady 


disc-wind like 


34 


640 


72(0.56) 


4.3(0.03) 


9.5(0.19) 


13 


steady 


disc-wind like 



(t) Self-consistently determined average temperatures at the outer boundary 
(t) Values in brackets are the standard deviations of the time series values. 



increase. The optical depth in polar directions (6 ^ 60°) are almost 
always less than 0.1. Even at larger polar angles (60° ^ 9 ^ 85°), 
r < 1. Only along a couple of 6 angles near the equator (~ 90°) 
where the relatively high density gas is prominent (c.f. Fig. |4}, r 
becomes greater than 1. From this analysis, we conclude that our 
assumption of optically thin UV radiation in our model is reason- 
able. Since the input radiation field in our model becomes smaller 
as 6 increases (c.f. equation O and it becomes very small along 
the equator, the effect of the small optical depths at larger 6 angles 
found here is minimal. On the other hand, if the angular distribution 
of the radiation field peaks in the equatorial direction, then the UV 
optical depths found here would make some effect on our results; 
however, this is not the case considered in this work. 

If, however, the optical depth becomes very large (r S> 1) in 
the polar region in which our input radiation field is strongest, we 
would need to cons ider the eff ect of attenuation and scattering of 
the photons (see e.g.lOstriker. M cKee. & Klein 1991; Murray e t al., 
ll994l ; IOhsuga et al.l2005h . In particular, the scattered photons from 
the polar region to equatorial region would redistribute the angular 
distribution of our input radiation field. Furthermore, the scattered 
photons would heat up the relatively low temperature accreting gas 
near the equator (c.f. Fig.|4l. This would especially change the scale 
height of the 'disc-like' accretion structure, and the corresponding 
mass-accretion rate. To include the effect of multiple scattering of 
photons, one would need a proper treatment of radiative transfer. 
Implementing such an effect in our model would require a major 
change in our computational code, and this is beyond the scope of 
this paper. 



5 CONCLUSIONS 

We have performed a set of 2-D (axisymmetric) hydrodynamical 
simulations (TableO to study radiation-driven outflows from a cen- 
tral part (< 10 pc) of AGN. We have extended the radiation-driven 
AGN outflow model of lProgal ( 120071) by relaxing the assumption of 
a constant accretion luminosity. This allows us to determine the ac- 
cretion luminosity consistently with the mass accretion rate at the 
inner boundary, and consequently the two quantities are coupled 
through the radiation field. Therefore, we made another step to- 
ward a more comprehensive self-consistent hydrodynamical mod- 
els with radiative feedback. Although the accretion luminosity is 
self-consistently determined in our model, the angular distribution 
of the radiation field is not self-consistently determined, but it is 
fixed while simulations proceed. The radiation always peaks in the 
polar direction and decreases significantly in the equatorial direc- 
tion. 



We have examined the dependency of mass outflow rates 
on the accretion luminosity (Section [3. Il l for three different outer 
boundary temperatures (To = 2 x lO", 2 x lO'^ and 2 x 10* K) 
by varying the outer boundary density (po). In the following, we 
summarise our main findings. 

We found a relatively strong correlation between the mass out- 
flow rates A/out at the outer boundary (r — Vo) and the Eddington 
ratio r (Fig.[T}. The outflows found in our models can be divided 
into three groups: 1. thermally-driven winds, 2. radiation-driven 
winds with a wide inflow, and 3. radiation-driven wind with a nar- 
row equatorial inflow (disc-wind like solutions). Out of 22 mod- 
els considered here, the majority of them (19 models) produces 
radiation-driven outflow/winds (Groups 2 and 3). Only 3 models 
produce thermally driven winds (Grope 2: Models 15, 16 and 17). 
The thermally driven winds are found only in the highest temper- 
ature models (To = 2 x 10* K). The disc-wind like outflows are 
found in 4 cases (Group 3: Models 6, 7, 8 and 9) which all have 
the lowest temperature at r = To (To = 2 x 10® K). No outflow 
is found for the models with F < 0.3 when the ambient tempera- 
ture is To = 2 x 10** K, and for the models with T < 0.2 when 
To = 2 x 10^ K. 

For the models with radiation-driven outflows, we found that 
the correlation between A/out and F can be described as a power 
law. The slope of the index of the power law (g) is relatively in- 
sensitive to the outer boundary temperature ro(Table|2j. Using the 
models from all three values of To (see above), the index of the 
power law is found as q — 2.0 (±0.1). This is very similar to 



those of the radiati on-driven stellar wind (e. g. Castor et"al]|l975h 
and disc wind (e.g. IProga et al .11 19981 ; |Progall l999h which predict 
q ^ 1.67. This confirms that a very similar power-law scaling of 
mass outflow rates with the luminosity of the continuum source 
works at extremely different mass regimes, namely A'h ~ 1 Mq 
(cataclysmic variables), ~ 10 M0 (massive stars) and ~ 10* Mq 
(AGN), which also have different types of spectral energy distribu- 
tions and the geometry of the radiation field. 

As in Papers I, II, and III, we find that the outflow is driven 
from an inflow. However, here we found cases of the inflow-outflow 
solution of an extreme form (Section 13. 3t for the lowest outer 
boundary temperature models (To = 2 x lO'' K) but with rela- 
tively high outer boundary densities (po 2 x 10"^'' gem"''). In 
these models, the inflow occur in a narrow zone near the equatorial 
plane which resembles a thin accretion disc, and the outflow oc- 
curs in a wide polar angles (0 < 9 < 85°). This type of the flow 
is very similar t o a wind driven by radiation from a luminous ac- 
cretio n disc (e.g. IProga et aljri99a; |Progj|l999l ; IProga & KallmarJ 
l2002h . In our models, the inflowing gas is very slowly rotating, 
so that the rotation velocity of the inflowing disc-like structure is 
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much smaller than the Keplerian velocity. The outflows arise from 
the disc-like structure near the equatorial plane in a wide range of 
the radii (1.4 x 10"^ pc < r < 7pc). 

For the models with a relatively high density set at r = To, 
a system reaches a steady state at F > 1 (see Table [T] and Fig.[T](. 
In other words, we find the models can still accrete matter even 
with a super-Eddington luminosity (up to F 4). This is a multi- 
dimensional effect. In our axi-symmetric model, the radiation field 
or the strength of the radiation force is not spherically symmet- 
ric but has the cos 6 dependency (equation I IQt. This allows gas to 
accrete near the equatorial plane even with a super-Eddington lu- 
minosity (F > 1), as demonstrated in Fig. [5] Based on this figure, 
the correlation between Afout and F (for the radiation driven out- 
flow cases) is expected to continue up to F ~ 10; however, we are 
not be able to confirm this due to numerical difficulties in handling 
very strong shocks that occur in the flows in the models with F > 5 
(for To = 2 X lO'^ and 2 x 10** K models). In the range of F val- 
ues we have explored, we do not observe a complete shutdown of 
the accretion flow even at a super-Eddington luminosity, which is 
expected in a spherically symmetric model. The accretion is still 
self-regulating, but not as a spherical model predicts. 

When we relax the assumption of a constant temperature at the 
outer boundary and self-consistently determine the temperatures 
there, the Afout-F relation becomes very similar to that of the mod- 
els with the lowest outer boundary temperature. To = 2 x 10® K 
(Fig. |7}- The high accretion rate models (F > 2) show disc 
wind like solutions, as seen in the high accretion rate models with 
To = 2 X 10** K (Section[33}. 

The following is a list of some steps for our model that shall be 
considered in the future: (i) include the effect of scattered photons 
which might enhance heating of the cold equatorial inflow, (ii) ad- 
just the continuum spectral energy distribution based on the mass 
and the luminosity of the system, and (iii) extend the computational 
domain up to a few 100 pc scale to capture the gas dynamics in the 
narrow line regions of AGN, and (iv) include dust. 

In addition, we shall incorporate the results fr om smaller scale 
simulations (e.g. disc-wind models bv |Progall99l) since not all the 
gas which accretes though our inner boundary (which is relatively 
large i.e. r = 0.014 pc) would reach the central SMBH due to the 
wind from the accretion disc. 
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